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Abstract 


The  effect  of  various  impurities  and  micro-alloying  additions  (B,  N,  C,  O,  Al,  Si,  S,  and  P) 
on  the  intrinsic  resistance  of  the  23  (1 1 1 )  grain  boundary  in  tungsten  has  been  investigated  using 
the  molecular  dynamics  simulation.  The  atomic  interactions  have  been  accounted  for  through 
the  use  of  Finnis-Sinclair  interatomic  potentials.  The  fracture  resistance  of  the  grain  boundary 
has  been  characterized  by  computing,  in  each  case,  the  ideal  work  of  grain  boundary  separation, 
the  mode  I  stress  intensity  factor,  and  the  Eshelby  ’  s  Fj  conservation  integral  at  die  onset  of  crack 
propagation.  The  results  obtained  suggest  that  pure  tungsten  is  relatively  resistant  to  grain 
boundary  decohesion  and  that  this  resistance  is  further  enhanced  by  the  presence  of  B,  C,  and  N. 
Elements  such  as  O,  Al,  and  Si,  however,  have  a  relatively  minor  effect  on  the  cohesion  strength 
of  the  23  (1 1 1)  grain  boundary.  In  sharp  contrast,  S  and  P  greatly  reduce  this  strength,  making 
the  tungsten  quite  brittle.  These  findings  have  been  correlated  with  the  effect  of  the  impurity 
atoms  on  material  evolution  at  the  crack  tip. 
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1.  Introduction 


Due  to  their  high  density  and  strength,  tungsten  (W)-based  alloys  are  widely  used  in  Army 
systems,  especially  for  antiarmor  applications.  However,  technically  pure  W  is  extremely  brittle 
and  its  ductile-^brittle  transition  temperature  (DBTT)  is  typically  as  high  as  300°  C.  To 
overcome  the  problems  associated  with  the  limited  toughness  ofW,  the  so-called  W  heavy  alloys 
(WHAs)  have  been  developed  [1].  Manufactured  by  liquid  phase  sintering,  these  alloys  are 
actually  metal  matrix  composites  consisting  of  hard  W  particles  embedded  in  a  relatively  mild 
nickel-iron  (Ni-Fe)  matrix.  Although  the  density  of  WHAs  is  only  slightly  lower  than  that  of 
pure  W,  the  alloys’  antiarmor  performance  is  quite  inferior  relative  to  that  of  depleted  uranium,  a 
competing  antiarmor  material.  This  difference  is  related  to  the  fact  that,  while  a  depleted 
uranium  penetrator  develops  adiabatic  shear  bands  upon  impact,  which  give  rise  to  a 
“self-sharpening”  effect  [2],  the  WHA  penetrators  “mushroom,”  which  leads  to  an  effective 
increase  in  the  penetrator  diameter  and,  in  turn,  to  a  reduction  in  the  depth  of  penetration.  Owing 
to  the  well-known  environmental  problems  associated  with  using  depleted  uranium,  significant 
resources  are  currently  being  directed  toward  developing  WHAs  with  improved  adiabatic  shear 
band  behavior.  The  usual  approach  is  to  significantly  modify,  or  find  an  alternative  to,  the  Ni-Fe 
matrix. 

A  different  approach  has  recently  been  initiated  by  Krasko  [3],  who  suggested  microalloying 
as  a  way  of  obtaining  a  “ductile”  W.  Being  a  body-centered  cubic  (bcc)  metal,  such  W  is 
expected  to  acquire  the  necessary  adiabatic  shear  band  behavior  since  the  latter  has  already  been 
observed  in  high-strength  bcc  steels  in  conditions  of  “marginal”  ductility  [4]. 

The  reduced  cohesion  of  grain  boundaries  is  frequently  cited  as  a  major  factor  limiting 
ductility  and,  in  turn,  the  performance  and  reliability  of  high-strength  metallic  materials  [5,  6]. 
Intergranular  anbrittlement  in  metals  is  usually  caused  by  impurities  segregating  to  the  grain 
boundaries.  Impurities  present  in  bulk  concentrations  of  only  10-100  parts  per  million  (ppm) 
can  result  in  a  dramatic  decrease  in  the  mechanical  properties  (primarily  ductility,  fracture 
toughness,  and  fracture  strength)  of  metallic  materials  and,  thus,  pose  significant  processing  and 
application  problems.  This  detrimental  effect  of  mere  parts  per  million  amounts  of  impurities  is 
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consistent  with  the  fact  that  only  a  few  parts  per  million  of  impurities  are  sufficient  to  saturate  all 
the  grain  boundaries  in  a  polycrystal  of  the  typical  grain  size  (10-100  pm)  [3].  The  sensitivity  of 
DBTT  to  the  grain  size  confirms  the  aforementioned  effect  of  grain  boundary  impurities.  That 
is,  the  smaller  the  grain  size,  the  larger  the  amount  of  impurities  required  to  saturate  the  grain 
boundaries  at  a  given  temperature;  hence,  at  the  same  bulk  concentration  of  impurities,  higher 
toughness  levels  and  the  lower  DBTT  are  expected  in  fine-grained  polycrystals. 

In  addition  to  the  smaller  grain  size,  increasing  the  purity  of  the  metallic  system,  in  general, 
increases  the  material’s  ductility  and  lowers  its  DBTT.  For  example,  the  DBTT  in  hi^-purity  W 
single  crystals  processed  by  electron-beam  zone  remelting,  involving  the  use  of  a  special 
impurity  gettering  procedure,  was  found  to  be  as  low  as  -196®  C.  ff  impurities  are  the  main 
cause  of  embrittlement  in  W,  gettering  the  impurities  is  the  obvious  way  of  ductilizing  this  metal. 
A  well-known,  though  extremely  expensive,  purity  enhanced  procedure  is  die  so-called  “rhenium 
effect”  [7,  8].  While  the  actual  mechanism  of  the  rhenium  effect  is  not  yet  completely 
understood,  it  is  believed  to  be  based  on  the  gettering  of  oxygen  [7,  8].  A  more  promising  way 
of  removing  “harmful”  impurities,  such  as  O,  N,  P,  etc.,  from  the  grain  boundary  is  by  gettering 
them  with  microalloying  additions  of  selected  elements,  e.g.,  B,  which  leads  to  the  formation  of 
compounds  such  as  boron  oxides,  boron  nitrides,  etc.  [9-11].  This  process,  however,  requires  a 
careful  control  of  the  material  chemistry  since  ductility  improvements  due  to  impurity  gettering 
can  be  expected  only  when  the  resulting  compound  precipitates  remain  fine.  Any  coagulation 
and  coarsening  of  the  precipitates  generally  results  in  an  adverse  embrittling  effect. 

Using  both  semiempirical  and  first-principles  calculations,  Krasko  [3]  carried  out  a 
theoretical  analysis  of  clean  W  grain  boundaries  and  grain  boundaries  containing  various 
impurities.  His  work  provided  some  important  results  regarding  the  effect  of  various  grain 
boundary  impurities  and  microalloying  additions.  For  instance,  impurities  such  as  H,  N,  O,  P,  S, 
and  Si  arc  found  to  weaken  the  intergranular  cohesion  in  W.  On  the  contrary,  the  presence  of  B 
and  C  was  found  to  enhance  the  bonding  across  the  grain  boundary,  thus ,  improving  the 
intergranular  cohesion.  Furthermore,  the  so-called  site-competition  effect,  where  the  species 
with  a  lower  energy  at  the  grain  boundary  tend  to  replace  the  species  that  have  a  higher  energy  at 
the  grain  boundary,  was  found  to  play  an  important  role  in  affecting  the  impurity  distribution  at 
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the  W  grain  boundaries  [3].  Among  the  atomic  species  analyzed,  B  was  found  to  have  the  lowest 
energy  at  the  grain  boundary  and  thus  would  tend  to  displace  other  impurity  atoms  from  the  grain 
boundary,  while  at  the  same  time  enhancing  the  intergranular  cohesion.  Based  on  these  findings, 
microalloying  of  technically  pure  W  with  10-50  ppm  of  B,  was  recommended  as  a  way  of 
enhancing  the  ductility  of  this  metal  [3].  Experimental  investigation  of  W  alloyed  with  B  in 
amounts  15  times  higher  than  the  one  recommended  by  Krasko  [3]  was  found  to  result  in  a 
significant  (150°  C)  drop  in  the  DBTT  [9-11].  Microstructural  analysis  further  revealed  the 
presence  of  relatively  coarse  boron  oxide  particles,  which  are  the  result  of  the  excessive  amount 
of  B  in  the  W.  This  suggests  that  the  addition  of  10-50  ppm  of  B,  which  is  sufficient  to 
completely  displace  oxygen  from  the  grain  boundaries  via  the  “site  competition”  effect  without 
an  excessive  formation  of  boron  oxide,  as  suggested  by  Krasko,  may  result  in  even  higher  levels 
of  ductility  and  fracture  toughness. 

In  the  present  work,  the  molecular  dynamics  method  was  used  to  explore  the  effect  of 
impurities  and  microaUoying  additions  (B,  C,  N,  O,  Al,  S,  Si,  P)  on  the  fracture  resistance  of  the 
E3  (1 1 1)  grain  boundary  in  W.  The  organization  of  this  report  is  as  follows.  Modifications  of 
the  “environment-sensitive  embedding”  (ESE)  energy  functions  for  various  grain  boundary 
impurities  in  W  are  discussed  below.  Generation  of  the  computational  crystal  containing  a  crack 
along  the  E3  (111)  grain  boundary  is  presented  in  the  section  “Computational  Crystal.” 
Computations  of  the  ideal  work  of  grain  boundary  separation  and  the  Fi  conservation  integral  are 
discussed  in  the  next  two  sections,  respectively.  The  section  “Results  and  Discussion”  contains 
the  results  obtained  in  the  present  study,  while  the  main  conclusions  are  drawn  in  the  final 
section. 

2.  Procedure 

2.1  Modification  of  the  Finnis-Sinclair  Functions  for  W.  Within  the  embedded  atom 
model  (EAM)  formalism  [12, 13],  the  total  potential  energy  of  the  system  is  given  as  the  sum  of 
two  enwgy  contributions,  the  interaction  of  each  atom  with  the  local  electron  density  associated 
with  the  remaining  atoms  in  the  system,  called  the  embedded  energy,  and  a  pair-wise  interaction 
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term  reflecting  the  electrostatic  interactions  between  the  atoms.  In  particular,  the  total  potential 
energy  is  written  as: 


«»  =  sj-Cp,)  +  I  Afe,).  (1) 

where  F,  is  the  embedded  energy  of  atom  i,  pi  is  the  electron  density  at  atom  i,  and  (Fy)  is  the 
pair-wise  interaction  between  atoms  i  and  j  separated  by  the  distance  Ry.  The  electron  density  at 
each  site  is  computed  from  the  superposition  of  spherically  averaged  atomic  electron  densities 
as: 

A  =  I.pjfeJ.  (2) 

pK*#)  in  equation  (2)  represents  the  atomic  electron  density  at  site  i  due  to  an  atom  at  sitey  at  a 
distance  Ry.  The  superscript  a  in  p“{r,j)  is  used  to  specify  the  species  of  the  atom  at  site  j. 
The  embedding  energies  and  pair-wise  interaction  functions  in  equation  (1)  are  generally 
determined  by  fitting  various  physical  quantities  of  the  material  such  as  the  sublimation  energy, 
the  equilibrium  lattice  parameters,  the  second-order  elastic  constants,  the  various  defect 
formation  energies,  and  the  zero-temperature  equation  of  state  [14]. 

The  Finnis-Sinclair  method  [15]  that  is  usually  called  “the  N-body  potential  approach”  has  a 
slightly  different  foundation,  but  the  total  potential  energy  can  stiU  be  expressed  in  the  form  of 
equation  (1).  The  first  term,  however,  has  a  different  meaning  since  it  originates  from  the 
tight-binding  theory  according  to  which: 

F{Pi)  =  -A^,  (3) 

where  A  is  a  species-dependent  parameter  whose  value  for  W  (1.896373  eV)  was  determined  by 
Finnis  and  Sinclair  [15]  and  p,  now  represents  an  environment-dependent  effective  atomic 
coordination  number,  or  an  effective  atomic  density  at  site  i.  Finnis  and  Sinclair  [15]  originally 
proposed  the  following  functional  form  for  p,: 
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and 

Pf  fe; )  =  0 ,  for  Ry  >R^,  (4) 

where  the  superscript  h  stands  for  “host  atoms.”  For  W,  Finnis  and  Sinclair  set  the  cutoff 
distance  Rcat  to  a  value  of  4.400224  A.  Finnis  and  Sinclair  [15]  further  proposed  the  following 
form  for  the  pair-wise  potential  function  0/,-: 

('^0 + +  <^2*9  ).  for  S  c , 
and 

P,;K)=0,  fOTRy>C,  (5) 

where,  for  W,  the  cutoff  distance  c  is  set  to  a  value  of  3.4  A,  and  the  quadratic  polynomial 
coefficients  are  cq  =  47.13465,  ci  =  -33.7665655,  and  C2  =  6.254199.  For  W  containing 
interstitial  impurities,  Krasko  [3]  introduced  a  modified  Fiimis-Sinclair  formulation  by  simply 
adding  a  term  to  equation  (1)  to  account  for  the  energy  of  impurity  atoms  in  the  host 
environment.  The  modified  £tot  is  then  written  as: 

=  '£F{pt)+  \'Zhjk,)+  (6) 

i  ^  iyj  k 

where  the  first  two  terms  on  tiie  right-hand  side  of  equation  (6)  are  given  by  equations  (3)-(5), 
while  the  third  term  represents  the  sum  of  the  ESE  energies  of  the  impurity  atoms  and  is  equal  to 
the  interaction  energy  between  the  impurity  and  host  atoms  and  is  the  electron  density  at  the 
site  of  impurity  atom  k  due  to  the  surrounding  host  atoms  and  is  given  as: 

pT  =  X  p?  =  Sp?fc,)'  w 

where  p?(R,y)  represents  the  electron  density  defined  by  equation  (4)  at  the  site  A:  of  an  impurity 
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atom  due  to  the  host  atom  at  site  j,  while  Rkj  is  the  distance  between  the  impurity  atom  at  site  k 
and  a  host  atom  at  site  7.  As  mentioned  previously,  the  ESE  functions  introduced  by  Krasko  [3] 
were  derived  in  the  spirit  of  EAM,  which  is  relative  to  the  free  atom  electron  density  at  the 
impurity  site  due  to  the  host  atoms  pf ,  while  the  host  atom  embedded  functions  are  expressed 

relative  to  the  effective  atomic  density,  pf.  To  avoid  possible  confusion  relative  to  the  meaning 
of  die  density  functions  used,  the  ESE  functions  were  rederived  in  this  report  relative  to  the 
effective  atomic  density  used  by  Finnis  and  Sinclair,  pf .  The  resulting  ESE  {p^ )  expressed  as 
fourth-order  polynomials  are  given  as: 

ESEip^)=  +  C,p^  +  C,  +  C3  {p^j\  (8) 

The  values  of  the  coefficients  Co,  Ci,  C2,  and  C3  for  various  impurities  and  microaUoying 
additions  are  given  in  Table  1,  and  the  corresponding  ESE  vs.  functions  are  plotted  in 
Figure  1. 


Table  1.  ESE  Energy  vs.  Effective  Atomic  Density,  p.  Function  Coefficients  for  Various 
Impurities  in  W 


ESE{e\)  =  Vp(Co  +  C,p  +  cy  +  C,p>) 

Co 

c, 

C2 

C3 

B 

-0.8926 

0.0099 

-2.548  X  10'^ 

7.833  X  10’^ 

C 

-0.6177 

0.0125 

N 

-0.9423 

0.0157 

7.834  X  10’® 

0 

-0.6202 

0.0134 

5.044  X  10‘* 

A1 

-0.6259 

0.0244 

fc’TTi‘rnT‘iJ 

9.926  X  10'® 

Si 

-0.5201 

0.0233 

-8.199  X  10’^ 

8.131  X  10'* 

S 

-0.1674 

0.0224 

8.522  X  10‘^ 

P 

-0.3655 

0.0299 

- 1.292  X  10'^ 
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Figure  1.  ESE  Energy  vs.  Atomic  Electron  Density  Functions  for  Various  Grain-Boundary 
Impurities  and  Microalloying  Additions  in  W. 

2.2.  Computational  Crystal. 

22.1  Generation  of  the  Grain  Boundary.  Simulations  were  carried  out  using  a  cylindrical 
crystal  whose  crystallographic  orientation  is  shown  in  Figure  2.  The  dimensions  of  the  crystal 
are  given  in  terms  of  the  number  of  (TT2),  (1 1 1),  (110)  bcc  interplanar  spacings  along  the  three 
principal  axes.  The  initial  diameter  of  the  crystal  was  approximately  9  nm. 

To  generate  a  S3  (111)  grain  boundary,  the  atoms  in  the  upper  part  of  the  computational 
crystal  (x2  >  0)  were  rearranged  to  produce  a  configuration  that  is  a  mirror  image  of  the  lower 
part  of  the  crystal  across  the  grain  boundary  plane  (xz  =  0).  The  resulting  stracture  of  the 
computational  crystal  can  be  represented  by  succession  of  the  (1 1 1)  bcc  planes  as: 

cbacbacbAbcabcabc, 

where  the  grain  boundary  is  marked  by  A.  This  grain  boundary  is  of  a  tilt  type  and  is 
characterized  by  an  inverse  coincidence  lattice  number  S  =  3.  Krasko  [3]  showed  that,  due  to  the 
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.  X2—  [111]bcc 


Figure  2.  Schematic,  Size,  and  the  Crystallographic  Orientation  of  the  Computational 
Crystal  Used  in  the  Present  Work. 

associated  lowest  energy  level,  an  impurity  atom  such,  as  B,  C,  N,  etc.,  is  most  likely  to  occupy 
an  interstitial  position  in  the  center  of  the  trigonal  prism  formed  by  6  W  atoms  in  the  Z3  (111) 
grain  boundary.  Consequently,  all  the  simulations  of  the  Z3  (111)  grain  boundary  containing 
impurities  were  done  under  the  condition  that  the  impurity  atoms  occupy  the  interstitial  sites 
associated  with  the  trigonal  prisms  along  the  grain  boundary  plane.  Note  that  for  a  typical  grain 
size  range  of  10-100  jtim,  the  introduction  of  interstitials  into  every  grain  boundary  interstitial 
site  corresponds  to  the  impurity  concentration  range  of  80-20  ppm. 

2.2.2  Grain  Boundary  Crack.  To  generate  a  crack  along  the  E3  (111)  grain  boundary,  all 
the  atoms  in  the  computational  crystal  were  displaced  from  their  initial  positions  in  accordance 
with  the  plane  strain  linear  elastic  solution  developed  by  Sih  and  Liebowitz  [16]  for  the  crack 
displacements  in  an  amsotropic  single  crystal.  While,  in  general,  the  single-crystal  solution  of 
Sih  and  Uebowitz  [16]  is  not  valid  for  interfacial  cracks,  the  fact  that  the  two  crystals  in  Figure  2 
had  identical  corresponding  crystallographic  directions  in  the  three  principal  directions  justifies 
the  use  of  this  plane  strain  procedure.  Under  a  pure  tensile  load  applied  in  the  xt  direction  and 
for  plane  strain  condition  along  the  X3  direction,  the  components  of  the  displacement  along  the 
three  principal  axes  are  given  as  follows: 
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where  r  and  B  are  the  polar  coordinates,  i.e.,  the  radial  distance  and  the  angle  between  the  radial 
line  and  the  jci-axis  in  Figure  2.  Functions  pj  and  qj(j=l  ,2)  are  defined  as: 


Pi  = 

+ 

^12 

~  ^16^1 

P2  ~  ^11^2 

^12 

~  ^16^2 

_  ^12^1 

^22 

~  ^26^1 

- 

II 

d 

^22 

~  ^26^2 

(10) 


where  Ji  =  /iiCUi  +  and  S2  =  fi2CC2  +  are  two  noncomplex  conjugate  roots  of  the  following 
characteristic  equation: 

^11  f^J  ^^6  (^^12  ^66 )  ~  ^^26  f^j  ®  •  (i  1) 


The  remaining  two  solutions  of  equation  (11)  are  the  corresponding  complex  conjugates  of  the 
first  two,  i.e.,  fJa  =  fii  and  /i4  =  /r*2-  atj  coefficients  appearing  in  equations  (10)  and  (11)  are 
related  to  the  elastic  compliance  constants  of  the  material,  Sy  and  for  plane-strain  are  given  as: 
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^11^33  ~ 

^13 

ai2 

=  a^i 

_  S12S33  S13S23 

S33 

> 

^33 

^22^33  ” 

<;2 

*23 

=  aei 

_  ^16^33  ~  ^13^36 

S33 

> 

^16 

S33 

S«S33  - 

^36 

= 

_  ^26833  —  813835 

S33 

9 

^26 

^33 

(12) 


For  a  given  material,  the  elastic  compliance  constants  used  in  equation  (12)  depend  on  the 
crystallographic  orientation  of  the  computational  crystal  being  examined.  Parameter  Ki 
appearing  in  equation  (10)  is  the  mode  I  stress  intensity  factor  whose  critical  value  associated 
with  the  reversible  crack  extension,  the  Griffith  stress  intensity  factor,  is  given  as  [16]: 


^Gr  = 


where  A  is  defined  as: 


A  = 


^11^22 


Nl/2 


“22 

\^n  j 


2ai2  +  agg 


2a 
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^1/2 


(13) 


(14) 


and  this  parameter  also  in  the  relationship  between  the  /i-integral  and  the  stress  intensity  factor 
,  /j  =  AKi  by  Sih  and  Liebowitz  [16].  The  ^ysep  term  appearing  in  equation  (13)  represents 
the  ideal  work  of  grain  boundary  separation  and  is  equal  to  the  amount  by  which  the  work  per 
unit  area  of  the  crack  surface  done  by  the  external  loads  exceeds  the  change  in  the  elastic  strain 
energy.  At  the  onset  of  crack  growth,  2ysep  is  given  as: 


=  r, 


+  rf 


(15) 


where  superscripts  A  and  B  denote  the  two  crystals  being  jointed  along  the  grain  boundary,  is 
die  surface  energy,  and  'VJnt  the  grain  boundary  (interface)  energy. 
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2.3  Computation  of  the  Ideal  Work  of  Grain  Boundary  Separation,  lymt-  To  compute 
the  surface  energy  associated  with  the  (lll)bcc  crystallographic  plane,  the  crack  plane,  a 
rectangular  slab  with  the  same  orientation  as  the  computational  crystal  shown  in  Figure  2  and 
with  its  top  and  bottom  surfaces  parallel  to  this  plane,  is  created.  The  height  of  the  crystal  is 
chosen  to  be  large  enough  to  avoid  the  interaction  between  the  top  and  the  bottom  surfaces. 
With  the  periodic  boundary  conditions  applied  in  the  other  two  orthogonal  directions,  the  energy 
of  the  slab  is  minimized.  The  surface  energy  is  then  defined  as  the  excess  energy  per  unit 
surface  area  of  this  crystal,  relative  to  the  one  of  the  same  size  in  which  the  periodic  boundary 
conditions  are  applied  in  all  three  directions. 

The  previous  procedure  is  repeated  using  a  slab  containing  tire  E3  (111)  grain  boundary  and 
the  corresponding  grain  boundary  energy  calculated  as  the  excess  energy  of  the  slab  per  unit 
grain  boundary  area  relative  to  the  perfect  bcc  without  a  grain  boundary  subject  to  the  periodic 
boundary  condition  in  all  the  directions. 

The  same  procedure  for  calculating  the  surface  and  the  grain  boundary  energy  was  utilized 
for  the  cases  where  the  initial  grain  boundary  contains  one  of  the  impurities  or  microalloying 
additions.  The  results  of  the  previous  procedure  for  pure  W  and  W  containing  impurities  are 
summarized  in  Table  2. 

2.4  Molecular  Dynamics  Method.  The  evolution  of  the  material  surrounding  the  crack  tip 
is  studied  using  the  standard  molecular  dynamics  procedure  at  100  K.  The  molecular  dynamics 
calculations  involve  the  solution  of  the  Newton  equations  of  motion  for  a  system  of  interacting 
particles  (atoms)  in  order  to  determine  the  classical  particle  trajectories  and  velocities. 

At  the  beginning  of  each  simulation  mn,  a  crack  was  generated  by  displacing  all  atoms 
according  to  equation  (9)  from  their  initial  positions  corresponding  to  the  perfect  unrelaxed  bcc 
biCTystal  structure  shown  in  Figure  2.  All  the  molecular  dynamics  simulation  runs  were  done 
under  the  fixed  periodic  boundary  conditions  in  the  X3  =  [llOjbcc  direction  to  ensure  the  plane 
strain  condition  in  this  direction.  Also,  at  each  level  of  Aj,  fixed  displacement  boundary 
conditions  were  applied  onto  the  outer  siirface  of  the  computational  crystal. 


11 


Table  2.  (lll)bcc  Surface  Energy,  75^,  S3  (lll)bcc  Grain  Boundary  Energy,  7gb,  and  the 
Corresponding  Ideal  Work  of  Grain  Boundary  Decohesion,  2'^,  for  W 
Containing  Various  Impurities 


Impurity 

7suxf(eV/A2) 

7gb  (eV/A^) 

27nt(eV/A2) 

pure  W 

0.4285 

0.5823 

0.2747 

B 

0.3006 

0.0909 

0.6382 

C 

0.3007 

0.1216 

0.6076 

N 

0.2583 

0.0943 

0.5925 

0 

0.2135 

0.3290 

0.3130 

A1 

0.1746 

0.3263 

0.2768 

Si 

0.1925 

0.3074 

0.3136 

S 

0.0735 

0.3515 

0.1505 

P 

0.1100 

0.3275 

0.2110 

The  temperature  at  which  all  the  simulation  runs  were  carried  out,  100  K,  was  set  by  initially 
assigning  to  each  atom  a  random  velocity  according  to  he  Boltzmann  distribution.  During  the 
simulation  process,  the  temperature  was  maintained  constant  by  exponentially  relaxing,  at  each 
time  step  (2  fs),  the  average  squared  atomic  velocity  using  a  time  constant  of  0.1  ps.  Using  this 
procedure,  die  temperature  could  be  maintained  with  3%  of  its  target  value. 

2^  Computation  of  Fi  Int^ral.  To  quantify  the  effect  of  the  crack-tip  stress  relaxation 
processes  on  the  resistance  of  a  Z3  (111)  grain  boundary  to  separation,  Eshelby’s  conservation 
F  integral  [17]  was  calculated  for  both  the  pure  W  case  and  the  cases  of  W  containing  one  of  the 
impurities  or  microalloying  additions.  The  F  integral  provides  a  means  for  determination  of  the 
energy  release  rate  accompanying  the  crack  extension  in  cases  where  material  nonlinearity 
effects  cannot  be  neglected.  F  is  a  vector,  and  its  components  along  the  three  principal  axes  (the 
cracking  direction  xi,  the  crack  plane  normal  direction  X2,  and  the  crack  front  direction  X3) 
represent  die  three  force  components  acting  on  the  crack  tip.  For  the  computational  crystal  and 
the  crack  orientation  used  in  the  present  work  (Figure  2),  Fi  acts  to  propagate  the  crack  tip,  while 
F2  and  F3  are  not  physical.  The  components  of  F  can  be  defined  using  a  closed  contour  F,  which 
surrounds  the  crack  tip  as  follows  [17]: 
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(16) 


F,  =  j;[wr.-<T.;^ 


dSj, 


where  W  is  the  strain  energy  density,  uj^Xi)  is  the  ^-component  (k  =  1,2,3)  of  the  displacement  at 
a  point  represented  by  the  coordinates  Xi(i  =  1,2,3),  Oij  are  the  stress  components,  and  dSj  =  dS.n;, 
where  nj  is  they  component  (j  =  1,2,3)  of  the  unit  outward  normal  vector  to  the  contour  segment 
of  length  dS.  To  calculate  the  stresses,  the  procedure  proposed  by  Hoagland,  Daw,  and  Hirth 
[18]  was  used.  According  to  this  procedure,  the  strain  energy  density  associated  with  an  atom  at 
the  position  Xi  is  defined  as: 


W'(x,)  = 


E(x,)-E, 

_  9 


(17) 


where  Eixi)  is  the  potential  energy  of  the  atom  with  the  coordinates  Xi  (i  =  1,2,3),  Eo  is  the 
equilibrium  energy,  which  the  atom  would  have  in  a  perfect,  stress-free  bulk  crystal,  and  Qo  is 
the  equilibrium  atomic  volume.  The  energies  E(xi)  and  Eo  are  calculated  using  the  Finnis- 
Sinclair  interatomic  potential  and  an  energy  expression  for  the  potential  per  atom  analogous  to 
equation  (1).  Assuming  that  for  small  strain  increments  the  strain  energy  density  can  be 
approximated  as  being  quadratic  in  incremental  strains,  the  stresses  are  calculated  using  the 
following  central  difference  formula: 


= 


dW  ^  Wjxi  +5i)-  Wjxi 
dey  2AUy 


(18) 


whore  Sy  are  the  components  of  the  strain  tensor  and  5,  =  uijxj  is  a  small  perturbation  of  the 
atomic  position,  which  is  the  result  of  the  application  of  a  small  uniform  displacement  gradioit 
Euy  to  the  entire  crystal.  The  value  of  Awy  =  0.0005  was  found  to  be  optimal  in  the  present  work. 
To  compute  the  components  of  the  stress,  one  of  the  Euy  is  set  to  0.0005  at  a  time  (the  other  Awy 
are  kept  at  zero)  and  corresponding  nonzero  S,  computed  and  used  in  equation  (18). 
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From  the  theory  of  elasticity,  it  is  known  that  the  distortions  uy  =  dui/dxj  are  composed  of 
strain  (e,y)  and  rotational  (wy)  componaits,  i.e.,  uy  =  ey  +  wy.  A  strain  ellipsoid  that  relates  the 
extensional  strain  between  an  atom  and  its  a-th  neighbor:  e“  =  (r“/R“)  -  1,  where  r“  and  are 
respectively  the  atomic  distances  after  and  before  straining,  can  be  used  to  obtain  the  local  strain 
components  e.y  as  follows: 


6“  =  6,, 


(19) 


where  are  the  directional  cosines  (in  the  unstrained  lattice)  to  a  atom.  Equation  (19)  simply 

represents  a  transformation  of  the  strain  components  from  tihe  coordinate  system  given  in 
Figure  2  to  the  one  whose  xi  axis  is  along  the  line  coimecting  the  two  atoms  in  question,  i.e.,  the 
atom  at  whose  location  the  strain  components  are  being  evaluated  and  its  a-neighbor.  To  find 
the  strain  components  ey  at  the  position  of  an  atom,  a  least-squares  procedure  was  used  to 
minimize  the  following  sum  over  the  n  neighbors  of  each  atom: 


4  = 


(20) 


The  minimization  procedure  dLJdey  =  0,  yields  six  linear  algebraic  (normal)  equations,  which 
were  readily  solved  using  the  Gauss  elimination  method  [19].  The  basis  for  determination  of 
rigid  body  rotations,  wy,  at  a  particular  site  is  obtained  by  proceeding  in  a  manner  similar  to  the 
one  employed  for  the  strains.  The  relationship  between  displacements  of  the  nearest  neighbors 
from  the  given  atom,  wf ,  and  the  antisynunetric  matrix  of  the  local  rotations  suggests  that  the 
following  sum  should  be  minimized: 

Lr  =  s[(<  +  W,2X“  -  Wi3X“)"  +  (u“  -  W,2X“  W23X“)" 

+  (u“  -H  Wi3Xf  -  W23X“)'  j.  (21) 
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The  minimization  procedure  ^LJ^/^Wij  =  0  yields  three  normal  equations,  which  were  also  readily 
solved  using  the  Gauss  elimination  method.  The  evaluation  of  the  Fi  integral  given  by  equation 
(16)  was  done  numerically  in  the  present  work  using  the  trapezoidal  rule  [19]  along  a  circular 
contour  around  the  crack.  Owing  to  the  discrete  nature  of  the  crystal,  all  the  atoms  within  a 
distance  of  0.2  nm  from  the  contour  were  assumed  to  be  associated  with  flie  contour.  The 
magnitude  of  the  contour  segment  dS“  corresponding  to  an  atom  a  is  calculated  as  dS“  = 
Ri^a  + 1  -  0a  + 1)/2,  where  R  is  the  contour  radius  and  Oa  + 1  and  Oa  - 1  are,  respectively,  the  polar 
angles  for  the  atoms  in  the  contour  that  immediately  follow  and  proceed  atom  a  in  the 
counterclockwise  direction.  The  unit  vector  normal  n“  to  the  contour  segment  dS“,  where  i  andj 
are  the  unit  vectors  in  xi  and  X2  directions,  is  defined  as  n“  =  cos  6a  i  +  sin  0a/,  respectively  [20]. 

Note  that  the  strain  energy  density  W  and,  hence,  Fi  each  have  two  contributions:  one 
associated  with  the  strain  energy  density,  and  the  other  arising  from  the  work  of  grain  boundary 
separation.  The  former  contribution  to  Fi  and  Ji  integral  acts  to  extend,  while  the  latter,  lysep, 
acts  to  close  the  crack.  Hence,  the  following  relationship  holds: 

Fi  =  (22) 

3.  Results  and  Discussion 

3.1  Ideal  Work  of  Grain  Boundary  Decohesion.  The  results  of  our  calculation  of  the 
(111)  surface  energy,  E3  (111)  grain  boundary  energy,  and  of  the  corresponding  work  of  grain 
boundary  decohesion  are  listed  in  Table  2.  The  results  pertain  to  the  case  when,  after 
decohesion,  all  the  impurity  atoms  remain  on  one  of  the  crack  faces  and  the  other  crack  face  is 
clean  of  impurities.  While  the  magnitudes  of  the  work  of  separation  take  on  different  values 
when  the  impurities  are  assumed  to  divide  equally  between  the  two  crack  surfaces  (data  not 
shown),  the  general  findings  regarding  the  effect  of  various  grain  boundary  impurities  on  the 
grain  boundary  decohesion  resistance  remain.  The  results  shown  in  Table  2  indicate  that,  on  the 
one  hand,  elements  such  as  B,  C,  and  N  inarease  the  work  of  S3  (111)  grain  boundary 
decohesion.  On  the  other  hand,  elements  such  as  O,  Al,  and  Si  have  a  relatively  small  effect  on 
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this  work.  Lastly,  elements  such  as  P  and  S  act  as  strong  grain  boundary  embrittlers  and  reduce 
tins  work.  These  findings  are  fiirdier  correlated  with  the  grain  boundary  energy  data  shown  in 
Table  2.  The  grain  boundary  energy  data  show  that  B,  C,  and  N  reduce  the  grain  boundary 
energy  the  most,  while  P  and  S  reduce  this  energy  the  least.  This  suggests  that  B,  C,  and  N 
would,  via  the  site-competition  effect,  tend  to  displace  P  and  S  from  the  grain  boundary  and,  in 
turn,  enhance  the  grain  boundary  cohesion.  In  other  words,  the  deleterious  embrittling  effect  of 
impurities  such  as  P  and  S  in  W  can  be,  as  previously  suggested  by  Krasko  [3],  reduced  by  using 
microalloying  additions  of  microalloying  elements  such  as  B,  C,  and  N. 

3.2  Initial  Atomic  Configurations.  The  initial  atomic  configuration  of  the  S3  (111)  grain 
boundary  in  pure  W  (“clean”  grain  boundary)  containing  a  crack  at  the  applied  stress  intensity 
factor  Kapp  =  2.0  Kqi  =  2.7694  MPa  m*^^  and  the  corresponding  configurations  after  relaxation 
(energy  minimization)  are  shown  in  Figure  3(a)  and  Figure  3(b),  respectively.  The  energy 
minimization  was  carried  out  using  the  conjugate  gradient  method  [21].  The  structure  of  the 
grain  boundary  ahead  of  the  crack  after  relaxation  can  be  characterized  by  analyzing  the  spacing 
of  the  adjacent  (1 1 1)  planes  as  a  function  of  the  distance  from  the  grain  boundary  and  comparing 
it  with  the  bulk  value  in  pure  W,  =  0.9137  A.  As  shown  in  Table  3,  the  (1 1 1)  interplanar 
spacing  oscillates  with  the  distance  from  the  grain  boundary  with  the  amplitude  of  oscillations 
gradually  decreasing.  Consequently  by  the  10-12th  plane  away  from  the  grain  boundary  (not 
shown  in  Table  3),  the  oscillations  are  practically  damped  out.  The  results  shown  in  Table  3 
correspond  to  the  average  (111)  interplanar  spacings  at  a  distance  between  16.9  and  19.6  A  from 
the  crack  tip  in  the  positive  xi  direction.  It  should  be  noted  that  the  distance  between  the  second 
and  the  third  plane  has  decreased  significantly  relative  to  that  in  the  bulk  bcc  W.  This  behavior 
is  the  manifestation  of  the  crystal’s  tendency  to  reduce  the  free-volume  of  the  grain  boundary. 
Also  note  that  because  the  plane  of  the  crack  lies  between  the  grain  boundary  plane  and  the  first 
(111)  plane  below  the  grain  boundary  plane,  the  variation  of  the  (111)  interplanar  spacing  with 
the  distance  from  the  grain  boundary  is  different  for  the  (111)  planes  below  and  above  the  grain 
boundary  plane. 

The  initial  atomic  configuration  of  the  23  (1 1 1)  grain  boundary  with  a  oack  in  W  containing 
B  impurity  atoms  located  in  trigonal  interstitial  sites  and  the  corresponding  configuration  after 
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(a)  '  (b) 

Figure  3.  Initial  Atomic  Configuration  of  the  ^‘Clean”  S3  (111)  Grain  Boundary  in  W 
Containing  a  Crack:  (a)  Before  Relaxation,  and  (b)  After  Relaxation.  The 
Applied  Level  of  the  Stress  Intensify  Factor  Aapp  =  2.0  Kqt  =  2.7694  MPa 

Table  3.  The  Spacing,  in  A,  Between  Two  Adjacent  (111)  Planes  as  a  Function  of  the 
Distance  From  the  Grain  Boundary  in  Pure  W  and  in  W  Containing  B  and  P 
Grain  Boundary  Impurities 


Adjacent  (1 1 1) 
Planes 

Type  of  Grain  Boundary  I 

Clean 

WithB 

WithP 

1  Above  the  Grain  Boundary  Plane  I 

G.B.  Plane-Plane  1 

1.292 

1.334 

1.555 

Plane  1 -Plane  2 

0.609 

0.801 

0.906 

Plane  2-Plane  3 

1.009 

0.828 

0.734 

Plane  3-Plane  4 

0.984 

0.987 

0.977 

Plane  4-Plane  5 

0.824 

0.865 

0.906 

Below  the  Grain  Boundary  Plane  | 

G.B.  Plane-Plane  1 

1.252 

1.338 

1.570 

Plane  1 -Plane  2 

0.566 

0.808 

0.908 

Plane  2-Plane  3 

1.051 

0.830 

0.734 

Plane  3-Plane  4 

0.973 

0.987 

0.968 

Plane  4-Plane  5 

0.821 

0.863 

0.917 
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energy  minimization  are  shown  in  Figure  4(a)  and  Figure  4(b),  respectively.  The  relaxed 
configurations  for  the  same  type  of  grain  boundary  containing  other  grain  boundary  cohesion 
enhancers  such  as  C  and  N  are  quite  similar  to  the  one  shown  in  Figure  4(b)  and  are  not  shown 
here  for  brevity.  The  variation  of  the  (111)  inteiplanar  spacings  with  distance  from  the  grain 
boundary  for  the  relaxed  configuration  depicted  in  Figure  4(b)  is  listed  in  Table  3.  In 
comparison  with  die  “clean”  grain  boundary  beyond  the  two  planes  adjacent  to  the  grain 
boundary  plane,  the  magnitudes  of  the  oscillation  amplitude  of  the  (1 1 1)  interplanar  spacings  are 
significantiy  lower.  In  addition,  the  distance  of  the  two  nearest  (1 1 1)  and  the  two  next-nearest 
(111)  planes  to  the  grain  boundary  increase  somewhat  due  to  the  presence  of  interstitial 
impurities  of  B,  C,  or  N  relative  to  those  in  the  case  of  clean  grain  boundary  in  W,  Table  3. 


(a)  (b) 


Figure  4.  Initial  Atomic  Configuration  of  the  £3  (111)  Grain  Boundary  in  W  Containing  a 
Crack  With  B  Atoms  (Smaller  Balls)  Located  in  Trigonal  Interstitial  Grain 
Boundary  Sites:  (a)  Before  Relaxation,  and  (b)  After  Relaxation.  The  Applied 
Level  of  the  Stress  Intensity  Factor  Kapp  =  1.4  Kgt  =  2.9548  MPa  m^. 

The  initial  atomic  configuration  of  the  £3  (1 1 1)  grain  boundary  with  a  crack  in  W  containing 
phosphorus  impurity  atoms  located  in  trigonal  interstitial  sites  and  the  corresponding 
configuration  after  energy  minimization  are  shown  in  Figure  5(a)  and  Figure  5(b),  respectively. 
The  relaxed  configuration  for  the  same  type  of  grain  boundary  containing  other  grain  boundary 
embrittlers  such  as  S  is  quite  similar  to  the  one  shown  in  Figure  5(b)  and  is  not  shown  here  for 
brevity.  The  variation  of  the  (111)  interplanar  spacings  with  distance  from  the  grain  boundary 
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Figure  5.  Initial  Atomic  Configuration  of  the  £3  (111)  Grain  Boundary  in  W  Containing  a 
Crack  With  P  Atoms  (Smaller  Balls)  Located  in  Trigonal  Interstitial  Grain 
Boundary  Sites:  (a)  Before  Relaxation,  and  (b)  After  Relaxation.  The  Applied 
Level  of  the  Stress  Intensity  Factor  Aapp  =  1.2  Acr  =  1.4563  MPa  m^. 

for  die  relaxed  configuration  depicted  in  Figure  5(b)  is  listed  in  Table  3.  In  comparison  with  the 
“clean”  grain  boundary,  the  magnitudes  of  the  oscillation  amplitude  are  significantly  larger. 
Furthermore,  the  distance  of  the  two  nearest  and  the  two  next-nearest  (111)  planes  from  the  grain 
boundary  is  significantly  increased  relative  to  those  observed  in  the  cases  of  the  clean  grain 
boundary  and  the  grain  boundary  containing  B,  C,  or  N  impurities. 

3.3  Material  Evolution  at  the  Crack  Tip.  The  progress  of  material  evolution  at  the  crack 
tip  in  the  case  of  a  clean  £3  (111)  grain  boundary  at  the  applied  stress  intensity  factor  Aapp  = 
2.7694  MPa  m^^^  is  shown  in  Figure  6(a)  and  Figure  6(b).  This  level  of  the  applied  K\  is  the 
minimal  level  at  which  a  clear  evidence  of  crack  propagation  by  0.5  ps  of  the  molecular 
dynamics  simulation  time  can  be  obtained.  A  comparison  of  the  atomic  configuration  after  0.5 
ps  simulation  time.  Figure  6(a),  with  the  original  configuration.  Figure  3(a),  shows  a  significant 
rearrangement  of  the  atoms  in  the  region  surrounding  the  crack  tip.  This  rearrangement  leads  to 
a  significant  crack  blunting.  No  evidence  of  dislocation  emission  from  the  crack  tip  can  be 
found.  Atomic  rearrangement  with  a  small  additional  crack  blunting  continues  to  take  place  with 


(a)  (b) 


Figure  6.  Atomic  Configuration  of  the  i:3  (111)  “Clean”  Grain  Boundary  in  W  After  the 
Molecular  Dynamics  Simulation  Times:  (a)  0.5  ps,  and  (b)  5  ps.  The  Applied 
Stress  Intensity  Factor  Aapp  =  2.0  Kcr  =  2.7694  MPa  m^. 

the  simulation  time,  which  can  be  readily  established  by  analyzing  the  atomic  configuration  after 
5  ps,  Figure  6(b).  It  is  interesting  to  note  that  as  a  result  of  the  aforementioned  material 
evolution,  a  one-atom  thick,  one-lattice  parameter-wide  strip  of  the  (llO)bcc  plane  has  formed  at 
the  very  crack  tip.  Figure  6(b).  This  plane  is  normal  to  the  original  crack  plane,  (111),  and  a 
[100]  direction  in  this  plane  is  aligned  with  the  initial  crack  firont-direction  [TlO] .  The  fonnation 
of  this  strip  not  only  causes  additional  crack  blunting  but  also  changes  the  local  ideal  work  of 
grain  boundary  decohesion,  which  in  this  case  can  be  approximated  as  where  is 

the  surface  energy  of  the  (100)  plane.  Using  the  aforementioned  procedure  for  the  calculation  of 
the  surface  energy,  the  ideal  work  of  grain  boundary  decohesion  has  been  determined  as  2^^^^ 

=  0.8128  eV/A^  Since  this  value  is  higher  than  the  one  0.2747  eV/A^  reported  previously  in 

Table  2,  this  may  be  the  reason  why  crack  propagation  stops  once  the  (110)  strip  forms  at  the 
grain  boundary. 


The  material  evolution  in  the  region  surrounding  the  crack  tip  in  the  cases  of  the  XS  (1 1 1) 
grain  boundary  containing  interstitially  dissolved  O,  Al,  or  S  atoms  is  quite  similar  to  that  of  the 
clean  grain  boundary  and  will  not  be  discussed  any  further. 
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The  progress  of  material  evolution  at  the  crack  tip  in  the  case  of  S3  (111)  grain  boundary 
containing  interstitially  dissolved  B  atoms  at  the  applied  stress  intensity  factor  iTapp  =1.6  Kot  = 
3.3770  MPa  m^^^  is  shown  in  Figure  7(a)  and  Figure  7(b).  K\  =  0.8  Kqi  =  1.6885  MPa  m*^^  is  the 
minimal  Ki  value  at  which  a  noticeable  crack  advancement  could  be  found  by  0.5  ps  of  the 
molecular  dynamics  simulation  time.  A  comparison  of  the  atomic  configuration  after  0.5  ps 
simulation  time,  Figure  7(a),  with  the  corresponding  original  configuration.  Figure  4(a),  shows 
that  a  significant  material  evolution  takes  place  in  the  region  surrounding  the  crack  tip.  This 
evolution  gives  rise  to  a  major  crack  tip  blunting,  but  no  evidence  of  dislocation  emission  can  be 
found.  The  aforementioned  material  evolution  and  the  accompanying  crack-tip  blunting 
continues  with  the  simulation  time.  Figure  7(b).  It  is  also  important  to  note  tiiat  after  5  ps.  Figure 
7(b),  B  atoms  form  a  cluster  at  the  crack  tip.  This  cluster  is  expected  to  enhance  the  intergrain 
bonding  and,  in  turn,  the  ideal  work  of  work  of  grain  boundary  decohesion.  This  may  be  the 
reason  why  once  the  cluster  of  B  atoms  is  formed,  the  crack  propagation  ceases. 


(a)  (b) 


Figure  7.  Atomic  Configuration  of  the  23  (111)  ^^Clean”  Grain  Boundary  in  W  After  the 
Molecular  Dynamics  Simulation  Times:  (a)  0.5  ps,  and  (h)  5  ps.  The  Applied 
Stress  Intensity  Factor  Aapp  =  1.6  Acr  =  3  J770  MPa  m^. 


The  progress  of  material  evolution  at  the  crack  tip  in  the  case  of  the  23  (1 1 1)  grain  boundary 
containing  interstitially  dissolved  P  atoms  at  the  applied  stress  intensity  factor  ^Tapp  =1.2  Kqi 
=  1.4563  MPa  m^^^  is  shown  in  Figure  8(a)  and  Figure  8(b).  Ky  =  0.4  Kot  =  0.4854  MPa  m*^^  is 
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(a)  (b) 

Figure  8.  Atomic  Configuration  of  the  S3  (111)  Grain  Boundary  in  W  Containing 

Phosphorous  Atoms  (SmaUer  Balls)  After  the  Molecular  Dynamics  Simulation 
Times:  (a)  0.5  ps,  and  (b)  5  ps.  The  Applied  Stress  Intensity  Factor  Aapp 
=  1.2  Kgt  =  1.4563  MPa  m^. 


the  minimal  value  of  the  stress  intensity  factor  at  which  a  noticeable  crack  extension  could  be 
observed  by  0.5  ps  of  the  simulation  time.  A  comparison  of  the  atomic  configuration  after  0.5  ps 
of  the  simulation  time,  Figure  8(a),  with  the  corresponding  initial  configuration.  Figure  5(a), 
shows  that  a  significant  rearrangement  of  the  atoms  takes  place  in  the  region  of  the  grain 
boundary  ahead  of  the  crack.  This  evolution  causes  a  significant  increase  in  the  interatomic 
distances  in  the  grain  boundary  region.  This,  however,  does  not  cause  any  crack  tip  blunting, 
and  the  crack  tip  remains  as  sharp  as  in  the  original  configuration.  The  aforementioned  material 
evolution  continues  with  the  simulation  time,  and  after  5  ps.  Figure  8(b),  the  crack  tip  has 
advanced  by  approximately  17  A  relative  to  its  original  position.  If  the  simulation  time  is 
increased  beyond  5  ps,  the  crack  continues  to  advance  until  it  reaches  fixed  atoms  at  the  outer 
rim  of  the  computational  crystal. 

3.4  Mechanics  of  Grain  Boundary  Fracture.  To  quantify  the  effect  of  the  various  material 
evolution  processes  described  in  the  previous  section  on  the  intrinsic  resistance  of  the  E3  (111) 
grain  boundary  in  pure  W  and  in  W  containing  different  impurities  microalloying  additions,  the 
Eshelby  Fi  conservation  integral  was  calculated.  As  discussed  previously,  Fi  represents  the 
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force  acting  on  the  crack  tip  in  the  xi  direction.  The  time  evolution  of  this  force  whose  initial 
value  corresponds  to  the  minimal  level  for  the  applied  stress  intensity  factor  at  which  a 
noticeable  crack  advance  could  be  seen  by  0.5  ps  simulation  time  is  given  in  Table  4  and  plotted 
in  Figure  9. 


Table  4.  The  Stress  Intensity  Factor  for  a  Reversible  Extension  of  the  Grain  Boundary 
Crack,  Kgt  [Equation  (12)],  the  Critical  Stress  Intensity  Factor  for  Crack 
Extension  Obtained  During  Molecular  Dynamics  Simulations,  £^crit>  and  the 
Corresponding  Value  of  the  Eshelby’s  Fi  Conservation  Integral 


jS:Gr(MPam^^^) 

^:crit(MPam''2) 

F,  (eV/A^) 

Ops 

msm 

Pure  W 

1.3847 

1.6616 

1.0953 

0.8397 

0.8256 

B 

2.1106 

1.6885 

1.4431 

0.6313 

0.6095 

C 

2.0587 

1.8528 

1.4006 

0.6127 

0.5913 

N 

2.0337 

1.9320 

1.3836 

0.6053 

0.5841 

O 

1.4781 

1.4781 

1.0055 

0.4424 

0.4269 

A1 

1.3900 

0.5560 

0.9456 

0.4137 

0.3952 

Si 

1.4795 

0.5918 

1.0066 

0.3857 

0.3722 

S 

1.0249 

0.4100 

0.6973 

0.5051 

0.4944 

P 

1.2136 

0.4854 

0.8257 

0.6612 

0.5485 

Figure  9.  Time  Dependence  of  the  Fi  Integral. 
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Li  each  of  the  cases  anal3^ed,  a  positive  value  of  Fi  is  required  to  initiate  crack  propagation. 
It  appears  that  the  initial  values  of  Fi  fall  into  three  groups:  the  highest  values  (Fi  ~  1.4  ±  0.1 
eV/A^)  for  the  E3  (111)  grain  boundary  containing  interstitially  dissolved  B,  C,  or  N  atoms;  the 
intermediate  value  (Fi  «  1.0  +  0.1  eV/A^)  for  the  clean  S3  (111)  grain  boundary,  and  the  one 
containing  dissolved  O,  Al,  or  Si  atoms  and  the  lowest  value  (Fi  «  0.7  ±0.1  eV/A^)  for  the  grain 
boundary  containing  P  or  S  impurities.  These  hidings  are  fully  consistent  with  the  results  for 
the  ideal  work  of  grain  boundary  decohesion  listed  in  Table  2. 

Figure  9  further  shows  that,  as  the  simulation  time  increases,  the  aforementioned  material 
evolution  at  the  crack  tip  continues  to  take  place,  and,  as  a  result,  the  magnitude  of  the  crack  tip 
decreases.  This  decrease  between  0  and  1  ps  is  significant  in  aU  the  cases  except  for  the  ones 
involving  interstitially  dissolved  P  and  S  atoms  at  the  grain  boundary  where  the  decrease  in  Fi 
with  the  simulation  time  is  quite  small.  Past  1  ps  no  major  change  in  the  magnitude  of  Fi  is 
found.  In  other  words,  the  material  evolution  at  the  crack  tip  causes  Fi  to  decrease  by  only  a 
lower  positive  value.  In  the  cases  of  the  clean  23  (111)  grain  boundary  and  such  a  boundary 
containing  B,  C,  N,  O,  Al,  or  Si,  this  lower  level  of  Fi  is  not  sufficient  to  give  rise  to  any 
additional  crack  extension,  and  consequently,  the  crack  growth  ceases.  In  sharp  contrast,  when 
the  grain  boundary  contains  P  and  S  impurities,  the  crack  growth  continues  while  the  Fi  remains 
fairly  constant.  The  continued  crack  growth  in  the  case  of  the  grain  boundary  containing  P  and  S 
atoms  is  promoted  by  a  continued  rearrangement  of  the  grain  boundary  atoms,  which  causes  an 
increase  in  the  interatomic  distances.  In  other  words,  the  atomic  density  in  the  grain  boundary 
region  continues  to  decrease  and,  in  turn,  lowers  the  cohesion  strength  of  the  grain  boundary. 

In  summary,  the  atomic  simulation  results  obtained  in  the  present  work  confirm,  as  has  been 
obsCTved  experimentally,  that  the  firacture  toughness  of  W  is  affected  greatly  by  the  type  of 
impurities  segregated  to  its  grain  boundaries.  Note,  however,  that  in  the  present  work  only  cases 
of  clean  grain  boundaries  and  boundaries  fully  saturated  with  a  single  species  are  considered. 
When  two  or  more  impurity  species  are  present,  Krasko  [3]  showed  that  the  one  whose  energy  at 
the  grain  boundary  is  lower  will  preferentially  segregate  to  the  grain  boundary  and,  thus,  have 
the  dominant  effect  on  materials  fracture  toughness,  hi  the  few  simulations  in  which  partially 
contaminated  grain  boundaries  (i.e.,  the  grain  boundary  was  not  fully  saturated  with  impurities) 
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were  used,  we  observed  that  the  character  of  the  effect  (embrittling  or  cohesion  enhancing)  did 
not  change  but  its  magnitude  was  lower. 

Note  that  our  results  differ  somewhat  from  the  findings  of  Krasko  [3].  Specifically,  N  was 
identified  by  Krasko  [3]  as  a  grain  boundary  embrittler,  while  in  the  present  work,  this  element 
was  found  to  be  a  cohesion  enhancer.  Furthermore,  O,  Al,  and  Si  were  identified  by  Krasko  [3] 
as  grain  boundary  embrittlers  while  these  elements  were  found  to  be  neutral  in  the  present  work. 
As  far  as  H  is  concerned,  the  embedded  energy  functions  for  this  element  were  not  available  and 
its  effect  on  grain  boundary  cohesion  could  not  be  studied.  We  believe  that  the  differences 
between  the  results  obtained  by  ICrasko  [3]  and  our  results  stem  from  die  differences  in  the 
simulation  methods  used — ^Krasko  [3]  applied  molecular  statics  and  constrained  the  motion  of 
atoms  to  the  planes  parallel  to  the  grain  boundary.  This  constraint  was  lifted  in  the  present  work, 
and  the  simulations  were  carried  out  using  molecular  dynamics.  This  method  allows  the  atoms 
to  relax  into  the  configuration  associated  with  a  lower  energy  and  is  generally  considered  as  a 
more  appropriate  simulation  method  for  nonzero  temperatures. 

4.  Conclusions 


Based  on  the  results  obtained  in  the  present  work,  the  following  main  conclusions  can  be 
drawn: 

(1)  With  respect  to  their  effect  on  the  intrinsic  fracture  resistance  of  the  S3  (111)  grain 
boundary  in  W,  the  atomic  species  analyzed  can  be  divided  into  three  groups: 

(a)  hiterstitially  dissolved  B,  C,  and  N  are  strong  grain  boundary  cohesion 
strengtheners,  which  both  enhance  intergrain  bonding  and  give  rise  to  a  lower  level 
of  effective  loading  at  the  crack  tip  by  promoting  crack  tip  blunting.  There  is  some 
evidence  that  clustering  of  the  interstital  atoms  is,  at  least  partly,  responsible  for 
the  observed  fracture  resistance. 
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(b)  O,  Al,  and  Si  do  not  significantly  affect  either  the  ideal  work  of  grain  boundary 
decohesion  or  the  crack  blunting  processes  relative  to  those  observed  in  the  case  of 
the  clean  13  (1 1 1)  grain  boundary. 

(c)  P  and  S  are  strong  grain  boundary  embrittlers  that  both  significantly  lower  the 
ideal  work  of  grain  boundary  decohesion  and  impede  with  crack  tip  blunting 
processes. 

(2)  Regardless,  whether  the  intrinsic  fracture  of  the  E3  (111)  grain  boundary  is  measured 
using  the  ideal  work  of  grain  boundary  decohesion  or  by  the  critical  value  of  the  Fi 
integral,  removal  of  deleterious  •  impurities,  such  as  P  or  S,  by  either  material 
purifrcation  or  by  the  site  competition  effect  with  B,  C,  or  N  can  have  a  great  potential 
in  achieving  the  goal  of  producing  the  “ductile”  W. 
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